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Abstract. To test the stability of the Parisi solution near T = 0, we study the 
spectrum of the Hessian of the Sherrington-Kirkpatrick model near T = 0, whose 
eigenvalues are the masses of the bare propagators in the expansion around the 
mean-field solution. In the limit T ^ 1 two regions can be identified. In the 
first region, for x close to 0, where x is the Parisi replica symmetry breaking 
scheme parameter, the spectrum of the Hessian is not trivial and maintains the 
structure of the full replica symmetry breaking state found at higher temperatures. 
In the second region T<g;a;<lasT— ^0, the components of the Hessian 
become insensitive to changes of the overlaps and the bands typical of the full 
replica symmetry breaking state collapse. In this region only two eigenvalues are 
found: a null one and a positive one, ensuring stability for T <g; 1. In the limit 
T — > the width of the first region shrinks to zero and only the positive and 
null eigenvalues survive. As byproduct we enlighten the close analogy between 
the static Parisi replica symmetry breaking scheme and the multiple time-scales 
approach of dynamics, and compute the static susceptibility showing that it 
equals the static limit of the dynamic susceptibility computed via the modified 
fluctuation dissipation theorem. 
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1. Introduction: 

The physics of spin glasses is still an active field of research because the methods 
and techniques developed to analyze the static and dynamic properties have found 
application in a variety of others fields of the complex system world, such as neural 
networks or combinatorial optimization or glass physics. In the study of spin glasses 
a central role is played by the Sherrington-Kirkpatrick (SK) model [1], introduced in 
the middle of 70's, as a mean-field model for spin glasses. Despite its solution, known 
as the "Parisi solution" [1113111], was found 30 years ago, some aspect are still far from 
being completely understood. In this work we discuss the spectrum of the Hessian of 
the fluctuations for the Parisi solution in the limit of vanishing temperature, a still 
not fully explored problem. 

The Hessian spectrum plays a central role non only for the stability of the Parisi 
solution of the mean-field SK model, but also for the study of finite dimensional 
systems. Its eigenvalues are indeed the masses of the "bare" propagators in the loop 
expansion about the mean-field limit. Thus the knowledge of the Hessian spectrum of 
the SK model is a prerequisite for any theory obtained from a development about the 
mean-field limit. 

The stability of Parisi solution for the SK model near its critical temperature T^,, 
has been established long ago 5, i6j by exhibiting the eigenvalues of the Hessian matrix. 
In few words, one has a Replicon band whose lowest eigenvalues are zero modes, and a 
Longitudinal- Anomalous (LA) band, sitting at (Tc — T) , of positive eigenvalues (both 
with) a band width of order (Tc — T)^. The analysis was partially extended later [7] 
via the derivation of Ward-Takahashi identities, showing that the zero Replicon modes 
would remain null in the whole low temperature phase, and hence would not ruin the 
stability under loop corrections to the mean-field solution. 

Despite these efforts a complete analysis of the stability in the zero temperature 
limit is still missing. Near Tc one can take advantage of the vanishing of the order 
parameter for T = Tc and expand the free energy, a simplification clearly missing close 
to zero temperature, where the order parameter stays finite. Moreover the T = limit 
is highly non-trivial. All these make the derivation of "effective" approximations valid 
for T ^- a rather diflicult task 013]. 

In this work, anticipating the main results, we show that in the limit T ^ 1 the 
spectrum of the Hessian can be divided into two regions. A first region where the 
spectrum maintains a structure similar to that found close to Tc, and a second region 
where only two eigenvalues, one null and one positive, are found. In the limit T — > 
the width of the first region shrinks to zero, and only the second region survives. 

The outline of the paper is as follow. In Section [H we describe how the Hessian 
of fluctuations associated with the SK model is obtained. In Section [3] we discuss 
the properties of the Parisi solution in the low limit T <C 1 and how these affect the 
Hessian spectrum by considering three simple cases. In Section |4l we show how spins 
averages, and response functions, involving any number of spins can be computed 
within the Parisi Replica Symmetry Breaking scheme with a finite number R of 
replica symmetry breaking steps. In Sections [5] and [6] using the results of Section 
m we derive the Hessian spectrum in the T — > limit for both the Replicon and 
Longitudinal- Anomalous Sectors. Finally Section [71 contains some discussions and 
conclusions. The two Appendices contain details on the calculation of spin averages 
in the continuous R oo limit, [Appendix A[ and the T ^ 1 limit, [Appendix B[ 
Finally in [Appendix C[ for completeness we report the approach in terms of frozen 
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2. Free energy functional, fluctuations and propagator masses 

The model is defined by the Hamiltonian [TU] 

H = s^Sj (1) 

where = ±1 are N Ising spins located on a regular d-dimensional lattice and the 
symmetric bonds Jij , which couple nearest-neighbor spins only, are random quenched 
Gaussian variables of zero mean. The variance is properly normalized to ensure a well 
defined thermodynamic limit iV — >■ oo. To average over the disorder one introduces 
replicas. After standard manipulations the free-energy density functional / in the 
thermodynamic limit is written as a function of the symmetric n x n site dependent 
replica overlap matrix Qf'' as [llj : 



/ n n ^Q"' exp£{Qf } (2) 



(ab) i 



P (ab) 



+ E In Tr,. exp (/^^ ^ Qf s^s') (3) 



(ab) 

where Qp^ is the spatial Fourier transform of Qf^ with respect to the site index i and 
/3 = 1/T. The notation "(ab)" means that sum is over distinct ordered pairs a < 6 of 
replicas. 

Equations ^ and ([3]) are the starting point of the perturbative expansion around 
the mean-field theory. One then writes 

gab ^ gab ^ ^gab (4) 

where Q""^ is the mean-field order parameter, and expands C in powers of SQ'f' , 

£ = £(0)+£(l)+£(2)^...^ (5) 

The first term 



£(0) ^ ^ 



2 



(6) 



E(Q"')' + liiTr,» exp(/?2^Q''^s'^^ 

(ab) (ab) 

gives the free energy density / in the mean-field limit, and equals that of the SK 
model. The second term reads 

£(1) = E E [Q"' - 

i (ab) 

where 

^ Tr,,.°.^ exp(/3^E(a,)Q°''^°^'') 

Tr,. exp(/32^(^,)Q-''s-s&) " ^' 

The vanishing of yields the stationary condition that determines the mean-field 
value of the order parameter Q'^^ = (s°s''), and ensures that tadpoles do not show 
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up in the loop expansion. Below the critical temperature the phase of the SK 
model is characterized by a large, yet not extensive, number of degenerate locally 
stable states in which the system freezes. The symmetry under replica exchange is 
broken and the overlap matrix Q'^^ becomes a non-trivial function of replica indexes. 
In the Parisi parameterization [l^ the matrix Q°'^ for R steps of replica exchange 
symmetry breaking is divided into successive boxes of decreasing size pr, with po = n 
and PR+i = 1, and elements given bjjjl 

Q"^ = Qr, r = 0,...,R+l (9) 

where r = aCi b denotes the overlap between the replica a and b, and means that a 
and b belongs to the same box of size Pr, but to two distinct boxes of size Pr+i < Pr- 
The solution of the SK model is obtained by letting i? — > cxd. In this limit the matrix 
Q'^^ is described by a continuous non-decreasing function Q{x) parameterized by a 
variable x, which in the Parisi scheme is x € [0, 1] and measures the probability for a 
pair of replicas to have an overlap not larger than Q{x). 

The meaning of x depends on the parameterization used for the matrix Q'^^ . In 
the dynamical approach |13j x labels the relaxation time scale t^, so that Q{x) — 
{s{tx)s{0)). Here the angular brackets denotes time (and disorder) averaging. The 
smaller x the longer t^. All time scales diverges in the thermodynamic limit but 
tx' /tx ~^ oo if X > x' . To make contact with the static Parisi solution one takes 
x € [0, 1], with X = corresponding to the largest possible relaxation time and x — 1^ 
to the shortest one. With this assumption one recovers Q(0) = and Q{1^) — qc{T), 
the largest overlap. In both cases Q(l) = 1, since it gives the self or equal-time overlap. 
Other choices are possible, e.g., those used in [TH [TSl [HI [HI [TB] to tackle the T ^ 
limit. We stress however that different choices just give a different parameterization 
of the function Q{x), but do not change the physics, since this is given by the possible 
values q that the function Q{x) can take and by their probability distribution P{q). 
This property is called gauge invariance [13l [HI [M] . In what follows, unless explicitly 
stated, we take for x the Parisi parameterization. 

The quadratic term 

P (ab) 

+ Y J2 6Qf5M''^-''''^6Qf. (10) 

(a6),(cd) 

defines the "bare" propagators of the theory. This quadratic form in SQf' contains 
the Hessian matrix 

= SfJbyxcd) - [{s's's^s'^) - (11) 

of the SK model whose eigenvalues rule the stability of the mean-field solution, and 
give the masses of the "bare" propagators. Terms with higher powers of SQf'' in the 
expansion (O defines the interaction vertices of the theory. 

In the reminder of this paper we shall consider the eigenvalue spectrum of 
the Hessian matrix M"'^'^'^ of the SK model for the Parisi solution in the very low 
temperature limit T <C 1. 

X The equality Q"'' = (s'^s^) that follows from the stationarity condition is valid only for a ^ b. For 
consistency one defines Q"" = Qij+i = 1. 
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2.1. The Hessian M"'''''^'^: the Replicon and the Longitudinal- Anomalous Sectors 

With 4 replicas the Hessian is characterized by 3 overlaps. We can distinguish two 
different geometries: 

(i) The Longitudinal- Anomalous (LA) Sector. This is characterized by the two 
overlaps r — a n b and s = c H d and, ii r ^ s, the single cross-overlap t = 
max[a n c, aH d,b n c,b n d]. Then we denote the matrix element in the LA Sector as 

M"''''"^ = Mt''\ r,s = 0,l,...R; t = 0,l,...R. (12) 

Note that < = i? + 1 if a = c or a = d or 6 = c or 6 = d. 

(ii) The Replicon Sector. In this case aOb = cDd = r, and the geometry is characterized 
by the two cross-overlaps 

u = maxfa f] c.aCid] ^ . , . „, 

7 7 ^ 71 u,v > r + 1 (13) 

V = max[6 n c, 6 n d] ' " ^ ' 

For the Replicon Sector the matrix elements are denoted as 

^a6;cd^^r;r^ W,W>r-fl. (14) 

The element M^:^, however, contains contribution from both the Replicon and LA 
Sectors, and one has [20] 

K;; = r^;; + K:^'' + K"'-^^'' (15) 

where the first is the Replicon contribution while the others come from the LA Sector. 
The latter can be projected out by taking the double Replica Fourier Transform (RFT) 
on the cross-overlaps u, v: 

R+l R+l 

^u = T.T. p-p^ - - + • (16) 



The LA terms indeed cancel in this expression and one can replace M^;^ in the double 
RFT by rM^^. This in turns implies that the inverse double RFT of Mr-r yields the 
Replicon contribution nM^]^ and not M^iJJ. 



3. How things work near T = 0: simplest cases 

The equation for Q{x) is rather difficult to solve by analytical and/or numerical 
methods for T — >■ 0. The origin of this difficulty can be traced back to the fact 
that, as the temperature decreases towards T = 0, the probability of finding overlaps 
Q"'' sensibly smaller than qc{T) = 1 - aT^ -\- 0{T^), with a = 1.575 . . ., vanishes with 
T [211 [22]- There is however a finite probability Xc ^ 0.524... that Q"'' < qc{T). 
Consequence of this the order parameter function Q{x) in the Parisi parameterization 
develops for T <C 1 a boundary layer of thickness 5 ^ T close to a; = 0, as shown 
in FigUl From the Figure we see that for very small T the function Q{x) is slowly 
varying for 5 <^ x < Xc. However, in the boundary layer < x < (5, it undergoes an 
abrupt and rapid change. In the limit T — > the thickness (5 ^ T — > and the order 
parameter function becomes discontinuous at a; = 0. 

Uniform approximate solutions valid for T ^ 1 can be constructed by using the 
boundary layer theory, that is by studying the problem separately inside (inner region) 
and outside (outer region) the boundary layer [23]. One then introduces the notion 
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Figure 1. Shape of the order parameter function Q{x) for T ^ 1 in the Parisi 
parameterization. The horizontal arrow shows the extent of the boundary layer 
of thickness <5 ~ T as T ^ 0. 



of the inner and outer limit of the sohition. The outer Umit is obtained by choosing 
a fixed x outside the boundary layer, that is in 5 ^ a; < 1, and allowing T — > 0. 
Similarly the inner limit is obtained by taking T — > with x < S. This limit is 
conveniently expressed introducing an inner variable a, such as a = x/S, in terms of 
which the solution is slowly varying inside the boundary layer as T — > 0. The inner 
and outer solutions are then combined together by matching them in the intermediate 
limit X — > 0, x/S — >■ oo and T — >■ 0. The inner solution Q{a) is a smooth function of 
a for T -> varying between and — 1 [Ml [Ml HB [H], similar to Q{x) at finite 
temperature. In the rest of this paper we shall concentrate on outer solution since as 
T ^> it covers the overwhelming part of the interval [0, 1]. 

The behavior of Q{x) for T ^ 1 has strong consequences on other relevant 
quantities, such as, e.g, the four-spin correlation entering into the Hessian matrix. 
We shall make this more quantitative in the next Sections. Here the only feature we 
whish to retain is that in the outer region for T <^ x < Xc and T — >■ 0, the function 
Q{x) is driven closer and closer to Q{xc) ~ Qc{T) as T approaches zero. It can be 
shown f55] , see also [Appendix B[ that for T ^ x < Xc and T — > 

Q{x)^l~c{(3xr^+(^^-o}j T^+0{{l3x)-\T'), (17) 

where c = 0.4108. . . and a — limT^.o(l ^ 'l{^c))/T'^- We note that the breakpoint 
Xc depends on T. The dependence is however very weak for low temperatures |22j 
and the approximation Xc(T) ^ Xc = 0.524 ... is rather good for T Q. From this 
expression we see that the variation of Q{x) in the outer region is 



Q{xc) - Q{x) ^ JT 
Q{x) 



1 



(18) 



so that one can safely take the approximation Q{x) ^ Q{xc) = qc{T) as T — )■ 0, the 
error being 0{T^) at least. Going back to R steps of Replica Symmetry Breaking this 
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(a) 



(b) 




(c) 



Figure 2. Tree configuration for replieas a, b, c with aHb = r. 



approximation translates into 

Qr-QR = qc{T) = l-aT^ + 0{T^), T^O (19) 

for all r in the outer region, that is such that T <C x{Qr) — Pr and T — > 0, or, 
equivalently, for fixed r ^ and T — > 0. We shall make this insensitivity with respect 
to the overlaps r in the T 1 limit more precise in the next Sections. Here we 
just discuss the consequence of the insensitivity on the elements of the Hessian by 
considering some simple cases. 

Suppose the two pairs of replicas are equal: (a, b) = (c, d). In this case from eq. 

one constructs the simplest Hessian component: 



i-{Q 



ah\ 



that for the overlap a H 6 = r gives 
Insensitivity implies that for fixed r and T — > 



wc have 



0. 



(20) 

(21) 
(22) 



The next simple case is when only three replicas are different, in which case we have 



b^c. 



(23) 



Ultrametricity imposes that the three replicas a, 6, c with aHb = r can be only disposed 
as shown in Fig. [2l The LA geometries (a) and (6) lead for T and fixed r and s 
to 

r — QrQs Q_r(1 — Qr) 



rp2 



T2 



-a + o(r"), 



(24) 
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while the Rephcon geometry (c) yields 

Q s - Ql ^ Qr{i - Q r) 

We shall see below that insensitivity implies that AfJ''' — A/^''' ^ 0, and that all 
Replicon components vanish. Then from eq. (|15p and p5|) it follows 



MZi;s = -^^W^ ^ ^-a + 0{T'). (25) 



M^f^ = -a + 0{T^), T^O. (26) 
Similarly from ([15]) and (|22p one obtains 

2M^f^ - M^'^ = l-2a + 0{T^) (27) 
which combined with (j26p gives 

Mj^'-^ = -1 + 0(T2), T^O. (28) 

The general case with four different replicas cannot be reduced to simple forms 
and the expression of the four-spin averages is required. This will be derived the next 
Section. 



4. Spin Averages 



The evaluation of the Hessian components requires the computation of the four-spin 
averages {s'^ s"^ s"^) for a generic geometry of the four replicas. This can be done by 
introducing the generating function 

Z{b) = exp{nG„i(6)} = Tr,. exp f 1 ^ A^f, s\,^ + ^ sA (29) 

\ ab a / 

where Aaf,, equal to jS'^Qab with /3 = 1 /T for the SK model, is a generic nxn symmetric 
matrix with Parisi's block structure, 

Spin averages follow from differentiation 



(s s 



lim 



1 



>o Z{b) 

Introducing the "block indexes" a/t, 
a = (ao,ai, . . . .an), 



d _d_ 

dba dbb 

Qk = 0, , 



Z{b) 



Pk 
Pk+i 



(30) 
(31) 



=b„=0 



- 1 



(32) 



where pk, with n = po > pi > ■ ■ ■ > pr > pr+i = 1, are the block sizes, the generating 
function can be written as multiple integrals over independent Gaussian variablesl^ 



» R 
Z{b) = / Vnia) n exp G,? (6f + 

■J n 



(33) 



where 2?/{(a) is the short-hand notation for: 

R „ t-i 



R 

= u n 



t—O ao,...,at-i 
§ We use Greek letters for summed replica indexes 



Dz, 



ao,...,at-i 



(34) 
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and = Zao^...^at-i ^-^^ independent Gaussian random variables of zero mean and 
variance one: 



Dz : 



dz 



27r 



The function G[i{b) is the "free energy" of a single spin in a field b: 

expG_R.(6) = Trs exp(6s) = 2 cosh 6. 
and the frozen (random) field , given by 
i?, 

where AAj — \t — A(_i, keeps track of the contributions from the various blocks. 

Inserting the form (155)) of Z{b) into eq. (I5T|) . and noticing that differentiation 
with respect to ba can be replaced by differentiation with respect to b^, we obtain 



(35) 
(36) 

(37) 



/ a b 
(S S ■ 



= lim 



1 



lim 



Z(0) 
1 



_d d_ 

db^dbf 



VR{a)\{e^^GB{b'^) 



oZ(0) 



2?K(«)nexpGfl(&f) 



(38) 



For any given geometry of the replicas a,b,c . . . the integrals can now be performed 
recursively from scale R up to scale 0. To illustrate the procedure let us consider 

R 



VRia)Y[e^pGRib^). 



(39) 



The field b^ can be written as 

6f - v/AA^zf +^yAA;z^ 



R-l 



= v/AA^^f + &f-"'- (40) 
Then splitting out the z^-integrals, and recalling that z^ depends only upon indexes 
ao, . . . ,aji-i, one has: 

„ R-i „ 

Z(0) = / 2?fl_i(«) n / Dz^exp[pj, Gni^AX;,z!^ + b^-')} . (41) 

This structure suggests introducing quantities Gr{b) as 

exp{p,G,_i(6^-i)} = j Dzl e^p[prGr{^f^rzl + ¥-^)] (42) 
so that eq. (p9| can be written as 



» R-l 

2(0) = / I)i?-i(a) n exp{p^Gfl_i(6«-i)} 



(43) 
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that has the same form of (|39p provided i? — > R—1. The entire process can be iterated 
up to level and leads to 

Z(0) = exp{poG'_i(0)} = J Dzel^p{poGoi^/\iz)Y (44) 

In the limit po ^ ji ^ one recovers the usual expression [12] 

G_i(0) = J DzGoi^/x'oz). (45) 

Equation has an interesting "physical" interpretation. The quantity Gii{h^) 
is the free energy of a system of one spin, i.e. of size 1, in the replica space in presence 
of the frozen field b^, that is with all random (Gaussian) z"^ held fixed. To move one 
level up, i? — i? — 1, we have to unfreeze and integrate over z^, while keeping all other 
fields z* with t < R frozen. The fields z* with t < R give the effective action, under 
the form of a (random) field, of the spins on the spin s" with a O b — t < R. Then 
integration over the field z^ means that only the spins s° and such that aOb — R 
are summed in the trace. All others are kept frozen. Thus the quantity Gii-i{b^^^) 
can be seen as the free energy (density) of a system in the replica space of size pr in 
presence of an external field b^~^, which gives the interaction with the frozen spins, 
that is the frozen degrees of freedom. Extension to the successive ^''-integration is 
straightforward. The quantity Gr-iib"^^^) is obtained by integrating out in turn the 
random fields z* with t > r, while keeping all z* with t < r frozen. This means that 
the trace is restricted to spins and such that a O b — t > r. The contribution 
from the spins not included into the trace, and hence frozen, is taken into account by 
the frozen field b^~^. The quantity Gr-i{b^~^) is then the free energy (density) of 
a system of size Pr in the replica space in presence of the external field b^~'^, which 
accounts for the degrees of freedom still frozen at scale r. 

The free energy G_i(0) is part of the total free energy density of the system, 
see eqs. ([6|) and (f29|) . and thus it is itself an intensive quantity in the real space. 
This implies that PrGr-i are intensive quantities, and hence as n — ?► the pr become 
densities in the real space: < Pr < 1. The Pr give a measure of the density of the 
frozen degrees of freedom at scale r — 1 as measured from the overlap. Consider indeed 
the function 

R 

x{q) ^ n + '^{pr+i - Pr)0{q - Qr) (46) 

which equals the number of pairs of replicas with overlap Q""^ less or equal to q: 
^il) = Pr+i if Qr < q < Qr+i- The fuuctiou x{q) is not decreasing with q, thus 
Pr < Pr' if r < r' as n ^ 0. Indeed in going from level r to level r—1 the number 
of unfrozen degrees of freedom, that is the number of spins in the replica space over 
which the trace is done, increases, and hence the number of frozen degrees of freedom 
decreases, as signaled by the decrease of the value of the overlap. This picture is fully 
consistent with the dynamical formulation of CHS [37] in terms of time-scales and 
density of frozen/unfrozen degrees of freedom. 

We can now turn to the problem of calculating spin averages. This differs from 
that of 2{Q) by the presence of terms that depends on the fields 6^, cfr. eqs. (pS)) 
and (p9|) . The recursion relation (j42|) is the usual rule to compute the free energy 
when some frozen degrees of freedom become unfrozen, and hence must be summed 
up in the trace. In the specific case those frozen at scale r but unfrozen at scale r—1. 
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represented by the fields z'' . The presence of Pr instead of Pr+i in the integrand fohows 
because at scale r there are Pr/Pr+i disjoint systems in the replica space, all with the 
same free energy, that merge at scale r — 1. This suggests the following recursion 
relation for the calculation of spin averages. Let Fs{b^) be a generic function of the 
field at scale s. We then define the quantity Fs-i{b^~^) at scale s — 1 as the average 
of -Fs(&a) O'^^'" the random field z" weighted with the statistical weight of the state, 
that is, 

i^,_i(6ri)=Z,_i(a)-i 



where 



X J DziZs{ar/P^+^Fs{VAX:zi + bi-^). (47) 

Zr{a) = Zr{bD = exp{pr+i Gr(6^)} 

= exp{p,+i GriV^Kzl + fe^-i)} (48) 
This recursion relation is supplemented by the boundary condition 

FM) =Fr{bl) (49) 

S—r 

where Fr{hr) is a known expression. Assume for example that 

Fr{Va) = ^Gr{bl) (50) 

then a simple calculation shows that 

Fr^iiK~')^T;^,Gr^r{V-^). (51) 
aba 

This result is not unexpected since it just states that the magnetization rn^ at any 
scale r is given by the derivative of the free energy of that scale with respect to the 
applied field at that scale: 

^rm^^Grm. (52) 

From this result it immediately follows that 

(0=m-iW = |^G_i(/.) 
_ d_ 

where h is an external field[jil Clearly (s°) — if h — 0. 

To compute the two-spin correlation {s°'s^) with a O b = r, i.e. the overlap Qr, 
we have to evaluate the integral 

/R 
Vnia) \{ ZR{a) mR{b^) mR{b^) (54) 
n 



DzGoiV>^oz + h)= DzmoWXaz + h). (53) 



where we used (j52p . On a branch-tree diagram the two replicas a and b with adb — r 
are on different branches for scales s > r, see Fig. [3l i.e., they belong to different 

II The field h includes a factor (5 from the statistical weight exp(— /SH), see also eq. Thus the 

correct expression in terms of the real external field would have fih. We prefer to leave the factor /3 
hidden into the field to have no factors /3 into the definition of the spin averages via eq. I)31|l . 
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b d 

Figure 4. Longitudinal- Anomalous Sector tree configuration. 

systems. The two local fields 6* and 6^, and hence the magnetizations, are independent 
from each other and the integration factorizes. Then 



I Vr{a)\{Zr{a) mr{blf (55) 

since at scale r the two replicas end up in the same system of size pr and the two fields 
and 6]J become equal. This expression suggests introducing the quantity 

Fs{K) = ^?}{K) (56) 
with the boundary condition 

m^^}{bl)=mr{hlf (57) 
in terms of which we have 



= mi^'L,{h)''^'^ Dzml^^„{^z + h) (58) 

ailb—r ' J 

('2) 

where m^Q is computed from the recursion relation (|47p with ([SG)) and ([57|) . 

For higher order spin correlations we proceed in a similar way. Consider for 
example the four-spin correlation 

R. 

{s-s's'^s^r,^' -> / VR{a) II Zfl(a) mfl(&f ) mfl(&f ) m«(&f ) mfl(6f ) (59) 



where the replicas a, b, c, d have the LA Sector configuration shown in Fig. The 
four replicas are independent from scale R to i, where the replicas b and c end up in 
the same system in the replica space and the fields b\ and 6' become equal. Thus 

{s'^s's^sX'" ^ f Vt{a)f[Zt{a) mt{bl)mt{blfmt{b'a). (60) 
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Moving up along the tree the surviving three rephcas remain in different systems up 
to scale r, where the replicas a and b eventually find themselves into the same system. 

(2) 

Then, by using the quantity rrirj introduced for the two-spin correlation, we can write 



/ Vr{a)f[Zr{a) mr{b:)m['^ib:)mrib:,) (61) 

<J n 



where m[^}(b^) is obtained from the recursion relation (|47p with the initial condition 

mp(^(6*) = mt{b*)'^. For the next step we observe that replicas a and d remain into 
different subspaces up to scale s. Thus by introducing the quantity 

Fkib'j = m'i';^^{b^) (62) 
with the boundary condition, 

mi';}}{b:) = mr{b:)mi'}{b:) (63) 
we can move up along the tree up to scale s, and 

{s'^s's^srt' ^ / V,ia)f[Z,{a) mi';^} (bl) mM)- (64) 
•' 

The last step from scale s to scale is now straightforward. We introduce the quantity 
Fkm^ml';l;l^ib';) (65) 
with the boundary condition, 

^tls']m)^m[';']{bi)mM) (66) 
so that the final result reads 

{s'^s''s^srt^ = n^;^:%ih) J Dzm^^;^^^^i^oz + h). (67) 

In a similar way, once the replica geometry is specified, one can compute spin 
averages involving any number of spins. 

The above results, valid for any finite R, are easily extended to the continuous 
case i? — ^ oo. In this case, since the values of Q""^, and hence those of A^, are bounded 
in a finite interval, the differences AA^ — as i? — > c» to account for an infinite 
number of values in a finite interval^ As a consequence the recursion relations are 
replaced by differential equations. In particular the recursion relation (I42p becomes 
the Parisi equation [12j, while eq. (|47p is replaced by the partial differential equation 

F{x,b) = ~^^ F"ix,b) + 2xmix,b)F'{x,b) , x<t (68) 
where m{x, b) = G'{x, b), with the initial condition 

F{x,b) ^ F{t,b) (69) 

x—t 

where F{t, b) is some known expression at scale t. As usual the "dot" and the 
"prime" denote partial derivative with respect to x and b, respectively. Details are in 
[Appendix A| 

Wc conclude this Section by noticing that this formalisms can be easily extended 
to calculate the (static) response of the system to external perturbations that act at 

5 One can allow for a finite number of "jumps" , that is points where does not vanishes as i? oo. 
One then gets mixed- type solutions as those found, e.g. in spherical p-spin models |30i i31| . 



Stability of the SK model near T — 14 

given scales. Let us denote by ej^ a small external perturbation acting only on scale r. 
We then define the static response as 



(70) 

£'■=0 

where (s")*^ is the spin average in presence of . The perturbation can be seen as 
an extra contribution to the frozen field b^. Then, since acts only on scale r, in the 
recursion relation we end up with Gr+i(-\/ ^-^r+i-^a^"^ + + *^a)- By expanding to the 
first order in ej^, we finally have 

{sY^j-DAa) n j Dz:^^ e^^{pr+iGr+i{¥+')] 

{ao,...,ar)^(ao,---,ar) 

X J i?z^iexp{p,+iG.+i(5^+i)+p.+ie>.+i(6ri)} m,+i(6^+i) (71) 

and a similar expression for the normalization Zgr-(p). By neglecting all unnecessary 
indexes, the recursion relation (|T7| is then replaced by 



mf{b) = exp G;' (&)| 

Dzexp |pr+iGr+i ( V AA^+i z + b) + Pr+i mr+i(v/AAr+i z + 6)| 

X m^+i(v/AA^+i z + 6) (72) 



with 

exp|pr+iG^^ (6)| = J Dz exp^^pr+iGr+ii\/ A\r+i z + b) 

+ pr+i e'^ "1^+1 (V z + 6) I . (73) 
Taking the derivative with respect to e*", and setting e*" — 0, leads to 

— mf(6) _^p.+i(mf^^i,,(6)-m.(5)2). (74) 
Then from eqs. ([701 and ([7T|) we hav43 

= Pr+l (Qr+l - Qr) ■ (75) 

The derivative of the normalization factor Zgr(O) gives a term proportional to {s"")"^, 
which vanishes in absence of external field. The expression (I75|) is the static limit of the 
modified Fluctuation Dissipation Theorem introduced by CHS [57] in dynamics. 



5. Replicon Sector 

The Hessian is a x symmetric matrix that after block-diagonalization 

[211 [33] becomes a string of (_R + 1) x (i? + 1) blocks along the diagonal for the LA 
Sector, followed by 1 x 1 fully diagonalized blocks, for the Replicon Sector. The 
diagonal elements in the Replicon Sector are given by: 

AR,(r;fcf) = M]:j. (76) 

To evaluate M^'^ we need the matrix elements M^iJ, eq. (jlip . that is the four-spin 
average {s"" s'^ s'^) for the Replicon Geometry shown in Figure [5] From the results of 

+ The perturbation e"^ includes a factor /3. This removes the factor /3 on the r.h.s. 
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Figure 5. Replicon Sector tree configuration. 



Section S] and [Appendix A| this is given, in absence of an applied external field and 
for R CO, by 



,(2,2) 



ml2.2)(o,o). 



(77) 



The function mu,v {x,b) is solution of a chain of partial differential equation of the 
form 



F{x,h) = -'- 



F"{x,b) + 2xm{x,b) F'{x,b) 



(78) 



where m(x,b) is the local magnetization at scale x in presence of the field b. 

In our case, starting from the bottom of the tree in Fig. [SJ we first have to solve 
eq. (|78| for F{x, b) = m{x, b) and boundary condition 



Fix,b) 



m(l, b) = tanh6. 



(79) 



The range of x is u < a: < 1 for the left branch of the tree, and v < x < 1 for 
the right branch. To proceed towards scale we have to solve next eq. ([75]) for 



F{x,b) 



,(2) 



[x, b) and initial condition 



F{x,b) 



_ (2) 



ml''>{t,b) = m{t,by 



(80) 



where t = u,v depending upon we are on the left or on the right branch of the tree. 
The range of x is either r < x < u, left branch, or r < a; < w, right branch, see Fig. [5] 
To accomplish the last step, < x < r, we finally solve equation (|78p for 
F{x,b) — mu,v {x,b) and initial condition 



F{x,b) 



= m52)(r,6) = mi^^(r,6)ml,^^(r,6) 



(2)f 



(81) 



While these equations are valid for any temperature T, we are interested into their 
solution in the limit T — > 0. Following Pankov |251 one can show, see [Appendix B[ 
that in the outer region x ^ T and T — >■ the solution of the partial differential 
equation ([78|) looses its explicit dependence on the scale variable a:. As a consequence 
for i? — > oo the matrix element MJ:^ becomes independent of u and v for all r such 
that pr — x{Qr) ^ T as T ^ 0, and this in turns implies that M^'J is independent of 
k and I. Thus by exploiting this insensitivity we conclude that for all r in the outer 
region 



k;l r+l;r+l \ 



1 



7?.— >oo 



0. 



(82) 



The second equality follows from a Ward-Takahashi identity [J. In the inner region 
the Replicon spectrum maintains its complexity. However its relevance becomes less 
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and less important as T approaches zero, and vanishes in the limit T — > when the 
thickness of the boundary shrinks to zero. The Replicon spectrum, similarly to the 
order parameter function Q{x), becomes then discontinuous at x = 0. 

6. Longitudinal- Anomalous Sector 

The LA Sector corresponds to the 1) x 1) diagonal blocks along the diagonal 
labeled by the index fc = 0,...,i?+l. The matrix element in each block turns out to 
be 

laM^ = hir) S^l + \mI^' Si'\ r, . = 0, . . . , i? (83) 
where Aj,(r) is a shorthand for 

r M:'"^ k>r + l, 

k;r+l 

Afe(r) = <^ (84) 
M'y: ^ fc<r + l, 

I. r+l:r+l 

and Si'''^ = pi'"'' - pi%, fc = 0, 1, . . . , i? + 1, with 

M^''^ is the RFT of the matrix element AI^'^ with respect the cross-overlap t, that is, 

R+l 

Ml^^ = Y.Pt''''^(MP'-MtL\) (86) 

t=fe 

with , if r < s, 

Pt t<r 

p^"^ ^ { 2pt r <t<s (87) 
Apt r < s < t 

For finite temperature the matrix elements are different in each block, however for 
scales k in the outer region the RFT M^''^ and Aj,(r) become insensitive to the value 
of fc as T — 0. All correspondent blocks are then diagonalized through the single 
eigenvalue equation^ 

s=0 

where Sg = Ps — Ps+i- Li the outer region the eigenvectors satisfy /'' 7^ if 
T <Si x{Qr) < Xc as T — > 0, and zero otherwise. The eigenvalue equation then reduces 
becomes 

R 

Ala^ = TM35V^.r, r = r,...,R. (89) 

4 H-j-l — 2. 

s—r 

* The boundary term t = in the RFT is proportional to po = a-nd vanishes for n — > 0. The next 
term is proportional to pi M^'" , since Mg'" = 0, and vanishes as ij — > 00, so the only term which 
survives is k = R + I. 
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Figure 6. Possible tree configurations corresponding to the Longitudinal- 
Anomalous Sector with r < s and t > r. 



where r is the lower bound of the outer region: x{Q-) — x^SasT^O. The diagonal 
Replicon contribution vanishes for i? — > cx), as ensured by the Ward-Takahashi identity, 
and does not contribute. This equation has two distinct solutions. The first 



for Ef=F^s/' ^ 0- The last equality follows from eq. and S ~ (5 ~ T as T ^ 0. 

In the inner region, where the LA spectrum maintains the RSB structure, the solutions 
are smooth functions of the inner variable even for T — )■ 0. For T ^ the thickness of 
the boundary layer shrinks to zero, and the eigenvalues ([5(1)) and (PT]) cover the whole 
LA spectrum, with a discontinuity at a; = 0. 

7. Summary and Conclusions 

In this work we have presented the analysis of the very low temperature limit T ^ 1 
of the spectrum of the Hessian for the Parisi solution of the SK model. It has been 
long known that in this regime two distinct regions of the interval x G [0, 1] can be 
identified according to the variation of the order parameter function Q{x) with x. 
We have shown that this has strong consequences on the the structure of the Hessian 
spectrum. In the first region x < S ^ T, where Q{x) varies rapidly from Q(0) = up 
to Q{S) ^ <lc{T) ^ 1, the spectrum maintains the complex structure observed close 
to the critical temperature Tc for the full RSB state. We can call this region the 
RSB-like regime. In the second region, T <C x < Xc with Xc ~ 0.575 . . ., where Q(x) 



Ala = 
for Z]f=-^s/^ = 0, and 



(90) 




{x - - 13' {1 - q,.{T)) 

(a-l) + 0(T), r->0 



(91) 
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is slowly varying, the Hessian spectrum has a completely different structure. Here the 
components of the Hessian matrix become insensitive to changes of the overlaps and 
the bands observed in the replica symmetry breaking regime collapse. In this region 
only two distinct eigenvalues survive: a null one and the positive one. This ensures 
that the Parisi solution of the SK model then remains stable as the temperature goes 
to zero. Remarkable is the occurrence of zero modes in both the Replicon and the 
LA Sectors. Null eigenvalues arise from Replicon geometry, with Ward-Takahashi 
identities protecting them. Note, however, that the zero modes arise also from LA 
geometry, that is without protection of the Ward-Takahashi identities. 

We observe that for T <C 1 the order parameter function is almost constant for 
T <^ X < Xc, the variation being indeed of order [Q{xc) — Q{x)]/Q{x) = O [{T/x)^). 
Thus in this region we have a marginally stable (almost) replica symmetric solution, 
that becomes a genuine replica symmetric solution in the limit T — 0, with self- 
averaging trivially restored. It is worth to remind that the stability analysis of the 
replica symmetric solution also leads to two eigenvalues, one of which is zero to the 
lowest order in — T (and negative to higher order), and the other positive. 

In the limit T — the region where the RSB structure of the solution is found 
shrinks to zero, and only the RS part survives. This feature, in a sense, brings about 
some perfume of conciliation between aspects of Parisi mean-field approach and of 
the droplet approach [551 [IS]- We stress, however, that in order to identify a genuine 
droplet behavior, corrections to the mean-field have to be studied in more details. 

Concerning the multiplicity of the eigenvalues we observe that in each Sector, 
Replicon and LA, one has to separate the contribution from the RSB-like and the 
droplet-like regions. The former is proportional to the width S of the region. Therefore 
in the limit T — the contribution from the RSB-like region vanishes, and one has 
the usual Replicon and LA multiplicities for the droplet-like region. 

Finally, we have developed a method to compute spin averages in replica space 
involving any number of replicas, both for a finite number R of replica symmetry 
breaking steps and for the continuous limit R ^ oo. This generalizes some special 
cases known for the continuous limit i? — >■ oo, and to our knowledge it is new. Moreover 
it sheds light on the interpretation of the replica symmetry breaking method and its 
relation with the dynamical approach. For example we were able to compute the static 
susceptibility and show that it equals the static limit of the dynamic susceptibility 
computed via the modified Fluctuation Dissipation Theorem. 
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Appendix A. Spin Averages: continuous case 

The expressions derived in Section |4] are valid for any finite number R of replica 
symmetry breaking steps. Here we shall address the limit R ^ oo, where the replica 
symmetry breaking become continuous. If the values of are bounded in a finite 
interval, as is the case of the SK model, in the limit R oo they must be dense and 

AA^ ^ as i?^ oo, r = 0, l,...,i?. (A.l) 
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so that the integral in the recursion relations can then be evaluated by expanding in 
powers of AA^. We observe that a finite number of "jumps", that is values of r where 
AAr remains finite as i? — > oo, are possible. One then gets mixed continuous-discrete 
phase described by a piecewise order parameter function, as those observed, e.g., in 
spherical p spin-glass models [30l[31]. We shall not discuss this case here. 

Let us first consider the recursion relation for the free energy Gr ■ In the following 
we drop all unnecessary indexes. By expanding G'r(\/AAr z + b) in eq. (|42[) in powers 
of AXr and integrating over the Gaussian variable z, a straightforward algebra leads 
to 

[g;'(6)+p,(g;(6))'1 +o{Axi) (a.2) 



Gr-lib) ^ Gr{b) + 



where the "prime" denote differentiation with respect the argument, i.e. the field b: 

d 



GUb) = g-^Gr{b) 



(A.3) 



Next we observe that Gr{b) is a function oi pr, thus to extract the non-trivial part of 
eq. (|A.2p as AAr — > we have to specify what happens to Pr — Pr-i as i? — oo. 

Suppose Pr — Pr-i does not vanish in the limit i? — )■ oo. In this case eq. (|A.2p 
implies that 

Gr-i{b) ^ Gr{b) as AA^ ^ but p.;. 7^ (A.4) 

that is, Gr{b) does not depend on the scale. This is what happens, for example, in 
the SK model for Xc < r < 1, the so called "plateau". If p^ — Pr~i —J' as AA^ — J> 
then, by defining 

Aa; 



Pr = X Pr — Pr-1 

and changing the notation to Gr{b) — G{x,b), equation (jA.2[) leads for AA^ 
the Parisi equation [TH [35] 

X{x) 



G{x,b) = -'- 



G"ix,b)+x{G'ix,b)y 



where the "dot" denotes the (left) derivative with respect to the scale x, e.g. 



G(x,b) = -^G(x,b] = lim 

ox Ax^O 



Gjx, b) - Gjx - Ax, b) 
Ax ' 



(A.5) 

> to 

(A.6) 
(A.7) 



We note that with this definition the validity of the partial differential equation (|A.6p 
can be extended to include the case (jA.4[) since the derivative is always well defined. 
As a consequence the initial condition simply reads, see eq. (I36p . 

G(x = 1^,6) = ln(2 cosh fe) (A.8) 

while (|^5|) becomes 



G_i(6) = DzG{0,y^X(fi)z + b). 



(A.9) 

For the recursion relation (|47|) we follows a similar procedure, that is we expand the 
r.h.s. of (|47p in powers of AA^ and integrate over the Gaussian variable z. This leads 
to 

AXr 



Fr^l{b) = Fr{b) + 



F;'{b)+2prG'r{b)F;{b) 



As before if Pr 



' Pr-i does not vanish as AA^ 

Fr-lib) = Frib). 



OiAXl). (A.IO) 

0, the recursion relation reduces to 

(A.ll) 
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If, however, — Pr-i — Ax — )■ as AA^ — )■ we end up with the partial differential 
equation 



F(x,6) = -^\^F"{x,b) +2xm{x,b) F'{x,b) 



x<t (A.12) 



where m{x, b) = G'{x, b) is the magnetization at scale x in presence of a field b, with 
the initial condition 



F{x,b) ^F{t,b) 



(A.13) 



where F{t, b) is a known expression at scale t. 

As an example consider the two-spin correlation {s°'s^) with aO b ~ r. From eq. 
it readily follows that 



anb—r 



(A.14) 



(2) 

where mj- {x,b) is solution of the the partial differential equation 



mi^\x,b) 



Xix) 



ml^'>"{x,b) +2xm{x,b)ml^^'{x,b) (A. 15) 

(A.16) 



for < a; < r and initial condition 

ml^^x = r,b) = m{r, 6)^. 

To our knowledge this equation was first derived by Goltsev |37| . 

For the four-spin correlation {s°'s''s''s'^)l'^ of Fig. 2] we have a similar expression: 



where mp^\'^' (x, 6) is solution of the partial differential equation 



,(2,14)/ 



Xix) 



,(2,1,1)' 



t.r,s 



for < a; < s and initial condition 

_(2,1,1) 



f 2 1 1) 

{x, 6) -f 2 cc m{x, b) ra\^ ^ (x, b) 



The function m\ ^ (x, b) is itself solution of the partial differential equation 



(A.17) 

(A.18) 
(A.19) 



m'i^;^\x,b) = ^^\m['^;^^''ix,b) + 2xmix,b)m'i^;^^\x,b)\ (A.20) 



for s < X < r and initial condition 

ml ^ {x — r,b) — ml (r, b) m{r, b). (A. 21) 

Finally m.^^ (x, b) is solution of equation (jA.lsp for r < x < t and initial condition 
m,[^\x ^t,b)^ m{t,b)^. 

We conclude this Section by noticing that if we take F{x, b) = m{x, b), then the 
differential equation (IA.12[) becomes the known differential equation [11] . 



m(x, 6) = m"{x, b) + 2x m{x, b) m'(x, b) 



(A.22) 



with the initial condition 

m(l",6) tanh6 (A.23) 

derived by Sommers and Dupont for the local magnetization. Numerical solution of 
these equations can be found by using, e.g., the method described in Ref. p2] . 
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Pi^^y)^-^ F"{x,y)+2l3xm{x,y)F'{x,y) (B.l) 



Appendix B. The Pankov scaling regime 

To discuss the Pankov regime T — )■ and T ^ a; ^ 1 we first perform the change 
of variable b = (3y into the partial differential equation (jA.12|) to make explicit the 
temperature dependence. In the new variable the equation reads 

2/32 

where the "prime" now denotes differentiation with respect to y, and the local 
magnetization is m{x,y) = /3~^ G'{x,y), see footnote page[TT] Following Pankov [25] 
we assume that the dependence on the local fields y is via the combination z = Pxy, 
that is, 

F{x,y)^F{x,z), z = l3xy (B.2) 

and similarly for m{x,y). The differential equation (IB.1[) then becomes 

~ \{x) r~ ~ 
X F{x, z)= ^—^ [F"ix, z) + 2m(a;, z) z) 

^zF'{x,z). (B.3) 

Pankov has shown that in the outer region the "tilded" functions do not depend 
explicitly on the scale variable x: F{x,z) = F{z). All dependence on scale, field and 
temperature enters via the combination Pxy. Pankov called this the scaling regime. 
From eq. (|B.3p it is clear that the Pankov scaling regime is only possible iff 

x'^ X(x) c 

— c ~ constant => X(x) — const. ^ (B.4) 

2 x^ 

in which case the partial differential equation (|B.3[) reduces to the ordinary second 

order differential equation 

F"(z) = -\- + 2 m(z)j F'(z). (B.5) 
Lc J 

In the SK model X{x) = (3'^ Q{x), where Q{x) is the order parameter function, then 
from (jB.4p it follows that in the outer region Q{x) has the form 

c 

Q{x) = const. - ,^ ,„ ■ (B.6) 
Equation (HH) now follows by imposing Q{xc) = qc{T) = 1 - aT^ + 0{T'^). 

Appendix C. Descending the replica tree: the frozen field probability 
distribution functions 

In Section 2] we have shown how spin averages can be computed using a bottom-up 
approach, that is starting from level R at the bottom of the tree and climbing up 
towards level at the top of the tree. A top-down approach is also possible. 

To illustrate the procedure suppose we have to compute the following average, 

/ ^ _ JVAamlZric^) 9r{b:) 

where gr{ba) is a generic function of the frozen field b^ at scale r. For example 
9r{ba) ~ ^r{ba)^ for the two-spin correlation, see eq. ([55]). The average (|C.1|) can be 
rewritten in the simple form 

(g) = / dyPr{y)gr{y). (C.2) 
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by introducing the frozen field probability distribution function at scale r 

Following the procedure outlined in Section H] we can integrate the Gaussian variables 
in eq. (|C.3p . ending with 

where 

ffr-i(o; y) = exp 



2 



Pr[Gr{y)-Gr-l{b)]y (C.5) 



V27rAAr i 2AAr 
Equation (IC.4p has the same structure as eq. (|C.f I) , thus we can write 

Priy) = / dy'Pr-i{y') gr-iiy';y) (C.6) 



where Pr-i{y) is given by eq. (IC.SP with r — > r — f . Inserting the expression (jC.SP 
into (|C.6[) leads to the recursion relation 

Pr[y)= j Z?zexp{p, [G^(?/) -G,._i(yAA^z + y)]| 

X Pr-^i^fKK z + y) (C.7) 

that gives the frozen field distribution at scale r one it is known at scale r — 1. The 
initial condition is specified at level 0, at the top of the tree, and reads 

My)-^c.JJ^\. (C.8) 



as follows from eqs. (jC.4p and (jC.5p . 

In the limit i? — >■ oo (y) — >■ P{x, y) and the recursion relation ()C.7P is replaced 
by the partial differential equation 

Fix, y) = ^ [P"ix, y)-2x [m{x, y) P{x, y)] '] (C.9) 
where m{x,y) = G'{x,y), with the initial condition at a; = 0: 

Equations (jC.9P and (jC.lOP where first derived by Sommers and Dupont [M] using a 
variational approach. We note that taking gr{b) = mr{b)'^ from eq. (|C.2p one recovers 
in limit R oo the Sommers-Dupont expression 

Q(a;)= y dyP(x,2/)TO(a;,y)2. (C.ll) 

The approach in terms of frozen field distribution functions can be generalized to 
deal with averages of quantities that depend on more then one local field. Suppose 
for example that gr — > 9r{bai bl) with a^b. In this case eq. (IC.2P is replaced by 

(5)= \ dy\dy2Pr{vuyi)gr{y\-,yi) (C.12) 
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where -Pr (2/1,2/2) is the probabihty distribution function of the frozen fields yi and 
t/2 lying on two different branches of the tree at scale r. This satisfies the top-down 
recursion relation 

Pr{yi,y2)= DziDz2 Y[ exp [Griyk) - Gr^iiV AAr Zk + yk)] | 



fc=l,2 



X P^_i(v/AA^2;i +2/1, v/AA^22 + y2). (C.13) 

The initial condition is specified at the branching point s < r where the two branches 
meet, and reads 

Pr(2/i,2/2) -%i- 2/2)^(2/1). (C.14) 

r—s 

It is easy to verify that the frozen field distribution functions obey the sum-rule 

J dy2Priyi,y2) = Priyi), Vr. (C.15) 

In the continuous limit the recursion relation (jC.13l) is replaced by the partial 
differential equation 

A(x) 

k=l,2 



P(a;,2/i,2/2) = ^^y^ ^ |(9feP(x, j/i, 2/2) - 2 x 5^ [m(x, j/^) 9fcP(a;, yi, 2/2)] } (C.16) 



where dk = {d/dyk)- The generalization to frozen field distribution functions of any 
number of independent frozen fields is straightforward. 
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